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Mechanistic home range models are important tools in modeling animal dynamics in spatially- 
. complex environments. We introduce a class of stochastic models for animal movement in a habitat 

' of varying preference. Such models interpolate between spatially-implicit resource selection analysis 

'n— «' ' (RSA) and advection-diffusion models, possessing these two models as limiting cases. We find a 

closed-form solution for the steady-state (equilibrium) probability distribution u* using a factoriza- 
^..ji tion of the redistribution operator into symmetric and diagonal parts. How space use is controlled 

' by the preference function w then depends on the characteristic width of the redistribution kernel: 

when w changes rapidly compared to this width, u* oc w, whereas on global scales large compared 
to this width, w* oc w^. We analyse the behavior at discontinuities in w which occur at habitat type 
' boundaries. We simulate the dynamics of space use given two-dimensional prey-availability data 

I and explore the effect of the redistribution kernel width. Our factorization allows such numerical 

simulations to be done extremely fast; we expect this to aid the computationally-intensive task of 
model parameter fitting and inverse modeling. 

O: 
6 

X) 

O^i Due to the influences of habitats on the availability of food, shelter, mates and risk of predation, patterns of animal 
space use are strongly influenced by the spatial distribution of habitat types across landscapes. A widespread technique 
for analyzing such relationships between habitats and patterns of animal space use has been resource selection analysis 
(RSA) [li II 0], in which the intensity of space use is assumed to reflect an underlying resource selection function 
^f) ' giving an individual's preference for the habitat type found at that location. However, implicit in this approach is 
. an assumption that all habitats are equally accessible to an individual regardless of its current position: no account 
\l is made of the individual's finite movement speed or the spatial geometry of habitats on the landscape on which it 
moves 

More recently, an alternative framework has emerged in the form of mechanistic home range models [ill [l^ . In 
' contrast to RSA which is largely descriptive, such models yield spatially-explicit predictions for patterns of animal 
space use in the form of a probability density function (pdf), by modeling the process of individual movement. 
Mathematically, the fine-scale behavior of individuals is treated as a stochastic (Markov) process [1, [l^, [H, |T3|, 
specifying the probability of an animal at a given location moving to a subsequent location during a given time 
interval. From this description one can derive, in the limit of small time intervals, a continuous-time partial differential 
5^ , equation (PDE) for the evolution of the pdf. For example, a recent analysis of coyote home ranges in Yellowstone 
■ [i^l used a "prey availability plus conspecific avoidance" (PA-f CA) mechanistic home range model to account for the 
observed patterns of coyote home ranges within the park. In this model, individuals exhibited an avoidance response 
to encounters with foreign scent marks and a foraging response to prey availability (individuals decreased their mean 
step length in response to increasing small mammal abundance in different habitats.) In such mechanistic models, 
inferences about long-term space use are usually made by evolving the continuous-time PDE in order to converge to 
its steady-state; this can be computationally time-consuming [Til . 

In recent work [To| , we developed a mechanistic home range model which reconciles these two main approaches by 
combining the concept of a spatial preference function w{x) with a stochastic model of fine-scale movement behavior. 
At each time-step the movement of an individual is governed by its relative preference for the local habitat surrounding 
its current spatial position i.e. the preference function restricted to a region of size L centered on the individual's 
current location. The length scale L has two roles: it is the typical (jump) distance per time step, but also can be 
interpreted as the distance over which the animal is able to perceive differences in surrounding habitats. This new 
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a) RSA b) Proposed mechanistic model 




FIG. 1: Schematic comparing conventional resource selection analysis (RSA) against our proposed model for space use, in two 
dimensions. Shaded regions show areas of habitat in which the preference function w{x) is constant, a) RSA. The probability 
of finding an animal in each region is proportional to the product of preference (resource selection function) Wj and area of 
the region Aj. There is no account taken of the animal's current habitat type or location within the habitat, b) Proposed 
MERSA model for space use. The animal responds only to the preference function within a distance scale L from its current 
location (black dot). Each time step the future location (white dot) is chosen randomly from a localized distribution biased by 
the preference function. 

model, which we will call a MEchanistic RSA, or MERSA model, is compared against traditional RSA in Fig. [1] and 
detailed in Section |TT1 In [l^] we were able to show that, in one dimension (ID) with a spatially smooth habitat 
preference function, using the Kramers-Moyal expansion (0, and App. A and E of (llj ) one may derive an an 
advection-diffusion equation for the expected patterns of space use, with advection and diffusion coefficients related 
to the parameters of the underlying stochastic movement model. We then showed that the resulting steady-state pdf 
could be determined analytically, giving a intensity of space use proportional to the square of the preference function. 

In this paper we expand and generalize this result to the more biologically relevant case of individuals moving in two 
space dimensions across landscapes which may include discrete habitat types giving rise to a non-smooth preference 
function. Specifically, in Section Hill we solve for the steady-state pdf of our discrete-time stochastic model of fine-scale 
movement behavior directly and analytically, for arbitrary preference functions, without resorting to the conventional 
procedure of taking the continuous-time limit p^ . (Note that this relies on a particular algebraic feature of our 
model; for a general redistribution kernel such an analytic solution is not available.) By doing so, we are able rapidly 
to compute exact steady-state space use patterns for the model for a full range of possible length scale values L, 
rather than being confined to the limit of small L inherent in a PDE-based approach. We also gain an understanding 
of the numerically-observed pdf behavior at jumps in preference function, which had remained a mystery [lo| . A 
consequence of our solution is that the computationally-intenstive task of solving an inverse problem to fit multiple 
model parameters can become orders of magnitude faster. We discuss such numerical implications in Section IlII Al 
and give CPU timings for our numerical examples throughout. 

For illustrative purposes, we choose a translationally-invariant exponential distribution of jump lengths, a kernel 
which has proven useful for modeling coyote foraging movement [111] . However, our analytic solution also holds for a 
generalization of the model of to spatially-varying diffusion coefficient (spatially-dependent length scale L) such 
as occurs in modeling the prey-density dependent foraging rate of wolves and coyotes [a, [l^, [3 ■ In Section IIV Al we 
explore the transition from small L (where the model tends to the advection-diffusion equation derived in [10|), to 
large L (where the model becomes equivalent to RSA). Here 'large' and 'small' are in comparison to the typical spatial 
scale on which the preference function changes. In Section HVAI we explore this numerically both in ID, and with a 2D 
preference function derived from discontinuous real- world small mammal abundance data appropriate for coyotes. In 
Section FlV Bl we show rapid and efficient numerical simulation of the time evolution of the pdf. Furthermore in Section 
|V]we analyse the behavior at a sharp discontinuity in preference function, and show that on spatial scales larger than 
L this effect may be approximated by effective matching conditions in coupled advection-diffusion equations. Finally 
in Section IVll we discuss implications and draw conclusions. 
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II. MECHANISTIC SPATIALLY-EXPLICIT (MERSA) MODEL 

The time-dependent pdf of an animal we will represent by u{x,t), where x £ £7 is the location, fl C K.'' represents 
the habitat region of interest, and d is the dimensionality of space (usually 1 or 2). Thus, in a ID setting, u{x,t)dx 
is the probability that at time t a given individual is to be found in the interval [x, a; + da:]. (Note in 2D we use x 
rather than x to represent location vector.) Its normalization is 

u{x,t)dx — 1 for all t (1) 



n 



Our MERSA model is an example of a Markov process. At each time step of length r the current pdf at time t is 
acted on by a fixed linear operator to get the pdf at time t + t. This is expressed by the master equation 

u{x,t + T)= / k{x,x')u{x',t)dx' for a; e (2) 
Jn 

Given an initial pdf u(x, 0) — uo{x) for all a; G il, by iteration the pdf at arbitrarily large future times (multiples of r) 
may be computed. Here the redistribution operator kernel k{x,x') is defined as the conditional pdf of an individual 
animal's location x, a time interval r into the future, given that its current location is x' [2d |. This is an uncorrelated 
jump or 'kangaroo' process [14] (note the animal has no memory beyond the fixed time scale r). Since a redistribution 
kernel is a conditional pdf it is everywhere non-negative and ('columnwise') normalized by 

k{x, x')dx = 1 for all a;' G (3) 



n 

Clearly for a given ecological situation the choice of r determines the form of k{x,x'). For example, shorter r may 
demand a smaller kernel width L simply because animal speed is limited. By the length scale L we mean the typical 
size of |a; — a;'| for which k{x, x') is significant. We will not indicate explicitly the dependence on t of the form of k. 
The appropriate value of r depends on the application; it needs to be large enough that successive animal relocations 
can be approximately treated as uncorrelated. Real-world location data collection technology also can be a factor if 
fine-scale model fitting is to be done. For coyotes, a typical value of r is 10 minutes 

We consider a spatial preference (resource selection) function w G L^(fl) which controls relative preference for 
each location in the domain. We represent unbiased ('preference-free') diffusive animal movement with a symmetric 
redistribution kernel (f)(x,x'), that is, 

(/)(a;, x') = (j>{x', x) for all x, x' G fl. (4) 
The kernel (f) obeys the Markov normalization fS]); from this and symmetry it follows that 

(a:,a:')da;' = 1 for all x e fl (5) 

n 

which is the statement that a constant pdf is invariant under redistribution by cj). Since a uniform density gives no net 
probability mass flow, we say that (an operator with kernel) is advection-free. Our MERSA redistribution kernel k 
is this advection-free jump kernel biased by the preference function, in other words, 

, w{x)(j){x,x') 

k{x,x) = -— 6 

z[x') 

with normalization function (easily seen to be required to satisfy JSl)), 

z{x') := / w{x")(t,ix",x')dx". (7) 
Jn 

For ([6|) to be meaningful we must have z > everywhere; it is sufficient that w > everywhere for this to hold, which 
we will assume from now on. Note that since (f){x,x') may depend independently on x and x' (barring the symmetry 
constraint), it may represent a spatially- varying (and also anisotropic) diffusion coefficient. The MERSA model is 
thus more general than that of [10'], which was restricted to the translationally-invariant case 

(l3{x,x') = 4>{x-x'). (8) 

Here ^(•) is a function of relative displacement alone, which limits one to a spatially-invariant distribution of step- 
lengths in the underlying stochastic movement model. 
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A. Limiting cases of the model 



We now discuss two limiting forms of MERSA. Firstly, consider the case (j){x, x') = l/vo\{il) for all x, x' E fl. This 
corresponds to preference-free redistribution to a uniformly-random location in fi, without regard to current location. 
^ then becomes 

, , /N w(x) 

Hx,x')= \ „ . 9 

w(x")ax" 

Since this is independent of x' , within a single time step (and for all future time steps) the master equation reaches 
its steady-state pdf u* := limt^oo u(-, t) given by 

u*ix) = CM^), with normalization constant C^^ ^ f ^ix")dx" . (10) 

Jn 

This is formally equivalent to a conventional (time-invariant) RSA model, with linear dependence on preference. This 
assertion is illustrated when is divided into regions j — 1 • • • to each of area Aj and constant preference function 
Wj, for then (llOp assumes the more familiar RSA form Uj = ^j'u^j/(X]i=i -^kWk) where Uj is the probability of being 
in region j jlj, li,, i7j, see Fig. [1^. In the case of an infinite domain such as 17 = R'' (in which case w alone delineates 
the habitat) no constant normalizable (p exists; however, the above result may be reproduced by considering the limit 
in which a the kernel (p becomes much wider than all spatial scales of interest in the habitat (L oo). Then (j){x^ x') 
tends to a constant for x and x' within the habitat, and (jlOp again follows. Thus for both cases above we will call 
this the L ^ oo limit, and state that in this limit our MERSA model degenerates to RSA. 

Secondly, consider the L — > limit where k tends to a diag onal kernel. In order for time evolution to take place at 
all, we must also take the limit r ^ 0. It is well known |l3l. IT^ that the correct way to balance these two limits in 
order to reach a well-defined diffusion coefficient is to choose the variance of the kernel k to scale as r. In the ID case 
of smooth w, and a translation-invariant kernel ([5]), we have derived {10] that in this limit our MERSA model gives a 
Fokker-Planck PDE with known advection and diffusion coefhcients. From this we showed that the steady-state pdf 
is 

u*{x) = C2w'^{x), with constant C^'^ = w^{x")dx" , (11) 

that is, quadratic in preference function. (The integral is bounded since w is smooth and in L^{fl)). We note that 
in this Fokker-Planck limit, our model is equivalent to that of [3] with the 'potential' function U(x) = —2D\ogw{x) 
and constant diffusion d(x) — D. 

We will see in Section IIV Al how these differing L ^ oo and L Q steady-state limits are reached in practice. 

III. ANALYTIC FORMULA FOR STEADY-STATE PDF 

The condition that it* be a steady-state pdf is that it be invariant under the master equation ([2]), in other words, 

u*{x)^ I k{x,x')u*{x')dx' ioi-al\xen. (12) 
n 



Our main result is the following claim. 

Proposition 1 A steady-state pdf for the model redistribution kernel @) is given by 

u*{x) — Cw{x)z{x), with constant C^^ = I w{x')z{x')dx' , (13) 



Jn 

where the function z is defined by 

The proposition is proved by substituting © and into the right side of then noticing that z cancels, allowing 
the simplification 

k{x,x')Cw{x')z{x')dx' = C / w{x)(t){x,x')w{x')dx' 
n Jn 



Cw{x) / (l){x' , x)w{x')dx' 



in 

Cw{x)z{x), 

u*{x), (14) 
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verifying (fT2l) . Crucially, it is the assumption that 4> is symmetric that allows us to proceed from the first to second 
line. Notice that once standard assumptions about ergodicity are satisfied, the steady-state u* is unique {e.g. see 
Doeblin condition in Q p. 396; in our ecological context this is satisfied since we can assume that there is always some 
randomness to animal motion, i.e. the kernel 4> is always somewhat spreading at each spatial location). Since w is 
assumed to be everywhere positive, so is u*. 

We now derive some secondary results on the structure of K , giving intuition into the reason for existence of the 
simple formula (|13p . They may be skipped if no further mathematical insight is required. We switch to operator 
notation, expressing (fT^ as u* = Ku* where if is a Markov operator. Recall that a Markov operator is an integral 
operator with non- negative kernel obeying ([3]), which can be expressed K"^l = 1 where 1 is the constant function and 
-fT"^ the adjoint operator (with respect to uniform measure). Our model ^ expresses the factorization 

K^W<^Z'^, (15) 

where $ is the (Markov) operator defined by the integral kernel (j){x,x'), and W and Z are the diagonal operators 
which multiply by the functions w and z respectively. Thus the structure is an operator with constant invariant 
measure sandwiched between two diagonal operators. Note that z = since $ is symmetric. Then we have. 

Proposition 2 With the assumption w > c everywhere, for some c > 0, our model Markov operator il5\} 

1. satisfies detailed balance, that is, 

k{x,x')u*{x') — k{x' ,x)u*{x) for all X, x' fl, (16) 

2. is self-adjoint with respect to the measure 1/u*, and 

3. has all eigenvalues real. 



The proof is as follows. We define U* to be the operator multiplying by the function u*; it can be written U* = ZW . 
Using (|15p gives KU* — W^W , explicitly symmetric. In other words 

{KU*f ^ KU*, (17) 

which is equivalent to detailed balance (|16p . Now we use (•, •) to indicate the L^(0) (real) inner product with respect 
to uniform measure, and (-lOi/u* with respect to measure 1/u*. Using ((T7]) and the boundedness of 1/u* in the 
middle step we have 

{a,Kb),/^, = i^,KU*^)^iKU*^,^) = iKa,b),^^, for all a, b e (ft) (18) 
u u* u u* 

which proves part 2. Part 3 immediately follows by self-adjointness. We remark that the simple formula for u* 
can now be seen as a result of the need to symmetrize the factorization p5|) . 



A. Implications for fast numerical modeling 

We discuss briefly why the above result is important in practice. The analytic steady-state pdf formula is 
special because, for a general redistribution kernel no such formula exists. In particular, there is no formula for u* 
for commonly-used mechanistic kernels such as those expressed in polar coordinates with exponential radial jumps 



and a von Mises angular distribution (Ch. 3 of [ll|). (In that work all steady-state distributions had to be computed 
in the advection-diffusion limit, often in a numerically-intensive fashion; see App. G of [ll| and @). In fitting 
model parameters to real-world location data, a larg e number of such steady-state pdfs must be found as part of the 
parameter-optimization process {e.g. Sec. 4.3 of [ll|). 

Numerical solution of u* for any redistribution kernel requires discretization of the domain into N degrees of freedom. 
In 2D the N required for acceptable accuracy can be large [e.g. 10''). Solving for u* given a general kernel is then an 
eigenvector problem involving a (possibly dense) N-hy-N matrix discretization of that kernel. The iterative solution 
of such large eigenproblems can be slow especially when diffusion rates are small. In constrast, the formula in 
our MERSA model bypasses this and requires only the computational effort of the single matrix- vector multiplication 
required to evaluate the (discretized) integral ([7]). This is an acceleration by orders of magnitude. 

There is a further numerical advantage to a special case of the MERSA model. Namely if is translationally- 
invariant (diffusion coeffcient is spatially constant) then the action of the (discretized) convolution operator $ may 
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FIG. 2: Steady-state pdfs u* for exponential jump pdf of (|19|) for an increasing sequence of L values. The preference function 
w (thin solid line) is chosen to be smooth on the left side of the domain, and discontinuous and oscillatory on the right side. 



be computed in time 0(iVln A^) via the Fast Fourier Transform (FFT) [T^, which for large N is much faster than the 
0{N'^) dense matrix-vector multiply. This further speeds up computing u* via ([7]). Finally, the time-evolution u{x,t) 
can now be computed much more efficiently. A single time-step of the master equation Q may be performed with 
0{N hiN) effort by using (|15p in a split- operator method: division by z, followed by FFT application of <f>, followed 
by multiplication by w. 

We illustrate these numerical techniques and advantages below. All CPU timings are reported using a single core 
of a 2GHz Intel Core Duo processor running MATLAB 7.0 in GNU/Linux. 



IV. NUMERICAL RESULTS FOR AN EXPONENTIAL JUMP KERNEL 



In this section we illustrate the predictions of our MERSA model for steady-state and time evolution of space use for 
an idealized ID preference function, and for a 2D preference function that is based on spatially-complex, real-world 
measurements of prey abundance in different habitat types. Throughout we consider a translationally- invariant 
kernel, choosing an exponential kernel with width L, which takes the form in ID 



\P\/L 



2L 



where p := x — x', and in 2D 



^(P) = '/'(|P|) = 



1 e-IPl/-^ 
27ri \p\ 



(19) 



(20) 



where we remind the reader that p is a vector in the latter. This 2D kernel has been used as a model of animal 
movement, for instance describing the fine-scale movement of coyotes with remarkable accuracy [llj. Note that the 
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FIG. 3: Small-mammal abundance for Lamar Valley region of Yellowstone National Park collected by Crabtree et. al. (un- 
published); see Ch. 7 of [ij]. Units are the total biomass density of small mammals in kg/ha with darker colors indicating 
larger values. The x- and j/-axes are in kilometers. The highest abundances of small mammals are found in the mesic grassland 
habitats. 

2D kernel is radially symmetric, and when integrated over angle it gives a distribution of jump distances r := \p\ 
which is exponential, (1/L)e~''/^. 

A. Steady-state 

Here we will assume the domain has periodic boundary conditions. It has the interpretation that the piece of 
habitat f2 in question is surrounded by similar (repeating) habitat, often a reasonable assumption. Mathematically 
this condition is achieved in the case of the unit interval il = [0, 1) by replacing (|19p by a sum over a few nearby 
'image' kernels, 

1 ^ 

^^P^^^ E (21) 

m=-M 

where M is chosen such that (j) is periodic to some high accuracy {e.g. I0~^). An analogous 2D image sum is used in 
the 2D case where f2 is a rectangle. 

Note that our formalism can handle other boundary conditions. For instance we model a standard zero-flux 
boundary in Section flV Bl We have also checked that the periodic steady-state results are very similar to those for 
zero-flux, except when L is of order the size of the whole domain. (In this case for zero-flux u* is affected by the 
effective boundary discontinuity in w, a complication which we will not pursue here). 

1. One- dimensional case 

In Fig. [5] we show u* computed via (jl3p and ([7]) for the exponential kernel, for an ascending sequence of L values, in 
the unit interval. Our example preference function w has been chosen to exhibit a variety of lengthscales: for x < 0.6 
it is smooth, corresponding to a gradation in habitat preference, whereas for a; > 0.6 it is piecewise constant with 
discontinuous oscillations between the values 1 and 2 corresponding to isolated patches of more favorable habitat. 
The shortest length scale of w is 0.015, namely the size of the smallest constant patch near x = 0.76. We see that 
for very small L, u* accurately matches w'^ for all regions apart from those with the most rapid w variations. This 
matches the expectation in Section fllAI 

A gradual transition is seen in the sequence of Fig. [2] As L becomes larger than a given feature, u* in the vicinity 
of that feature starts to become locally proportional to w. Finally in plot (f), L is larger than any feature and u* 
globally becomes linear in w. The explanation for this transition is simple and lies with p3l) combined with the 
realisation that the function z is given by the function w smoothed locally over a width of about L. Fine (< L) 
features in w will thus be smoothed away giving a locally constant z, whereas for coarse (> L) features z « w and 
quadratic dependence in w results. This effect is common to any jump kernel with characteristic width L: fine-scale 
habitat features result in u* (x w, in accordance with the L — > oo limit pop . whereas coarse-scale habitat features are 
tracked according to u* oc w'^, in accordance with the Fokker-Planck limit (fTTI) . 
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Since u* is computed analytically via (|13p . its numerical accuracy is limited only by quadrature of the integral. We 
computed z using ([7]) via FFT convolution (Section IIII A[) in a few thousandths of a second on a uniform quadrature 
grid of = 400 points. This contrasts the order 1 s needed to iteratively solve for the dominant eigenvector of the 
dense matrix discretization of $ that would be required if no simplifying factorization of K or analytic solution for 
u* were known. 

2. Two-dimensional case with real-world habitat heterogeneity 

We start with the small-mammal density data B{x) shown in Fig. [3] where location x was sampled on a 0.1 km 
grid over a 7 km by 12 km domain Q. This density is piecewise constant, being derived from measured prey densities 
(mice, ground sq uirrels, pocket gophers and red-backed voles) appropriate for coyotes in six different habitat types 
(see Ch. 7 of [ll|). The main habitat feature is a strip of mesic grassland (the darkest region in the figure), which 
follows a valley floor and supports a high abundance of small mammals. 

We use a linear relationship between prey availability and preference 

w{x) ^1 + aB{x) for aU a; £ (22) 

where a (units of ha/kg) controls the strength of the preference per unit of prey biomass. The resulting steady-state 
pdfs u* are shown in Fig. [H This shows the patterns of space use for 3 different values of a, in combination with 
3 different values of length scale L ([20|) . Each column represents an a value, ranging from a weak preference (left 
column) to strong (right column). For comparison, the w function is shown at the top of each column. 

The computational grid was moderately-sized [N = 8591), the same spatial resolution as the underlying estimates 
of prey abundance was sampled. Once 4'{p) had been evaluated on the grid, solving for each steady-state pdf required 
only 0.009 s using with z computed from w by 2D FFT convolution 21]. This is 100-1000 times faster than 
an iterative solution for the dominant eigenvector of K if neither the factorization nor analytic formula are used (a 
single dense matrix application of a general K takes 0.45 s and many such iterations are required for convergence, 
the number depending on L and the particular w). Even if the factorization (split operator method) were used to 
perform each iteration, our analytic formula (I13p would still be 10-100 times faster. 

How is steady-state space use controlled by a and L? Comparing the L = 6 row (cases b,f,j) to the preference 
function itself (a, e, i), we see that with this large length scale u* is very close to proportional to u;, as in ID and as 
as explained in Section III Al Proceeding down the figure, we see smaller L values result in a u* with an exaggerated 
tendency towards space use becoming increasingly concentrated in areas of higher preference. This results in much 
more relative animal concentration in the mesic grassland region than would be predicted by traditional RSA. This 
tendency has reached its limit by the bottom L row (d, h, 1), where u* is close to proportional to lu^ . 

Consider the right-hand column (j, k, 1). Changing from L = 6 to L = 0.7 causes a substantial increase in relative 
space use in the western part of the mesic grassland (see the large bump to the left in k), but very little change in 
the eastern part of this same habitat type. The explanation is simple: the grassland strip is generally wider than 
L — 0.7 on the western side resulting in a local tendency towards u* oc u;^, but narrower than L — 0.7 on the eastern 
side giving here u* oc w. Finally in case 1 L = 0.1 is narrower than all parts of the mesic grassland and the densities 
equalize on west and east sides. This is analogous to the transition discussed in Section IIV A II This interesting 
geometric effect is absent in traditional RSA. 

B. Evolution in time 

In Figure [5] we show the time-dependent evolution of space use u{x,t) under the master equation ([2]) in a single 
space dimension for simple piecewise linear model preference function with a single discontinuity (see panel a), 
with exponential kernel with length scale L — 0.01. Zero-flux boundary conditions were created by using the non- 
periodic jump kernel (jl9p . Panels b-e show the resulting dynamics of space use, starting from an initial condition 
uo{x) = 6{x — xo) where xq ~ 0.82 

The computation with N — 400 grid points took 0.0005 s per time step; this was done with iterated multiplication 
by the dense K matrix since the non-periodic choice of cj) makes the $ operator no longer a convolution. Notice the 
discontinuity feature in w at a; = 0.75 establishes itself rapidly and persists throughout the full time range. 

Fig. [HI shows the time-depednent evolution of space use for the 2D real- world biomass preference model (|22p with 
a = 500 and L = 0.7 (same parameters as Fig. |3k), for a delta-function initial condition at xq = (9,4). After a 
single iteration (panel a), the local preference biases the individual's movements towards the narrow nearby eastern 
strip of mesic grassland. In panel (b) its space use is still split between an expanding radial distribution about its 
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initial position and the nearby grassland. In panel (c), the compounded mechanistic movement steps have caused the 
individual's space use to become concentrated in the eastern strip of mesic grassland but the mesic grasslands in the 
western portion of the landscape are, at this stage, mostly unoccupied. Much later in the simulation however (panel 
d), the intensity of space use in the western mesic grassland areas is higher than the eastern mesic grasslands, due to 
the geometric effects on the steady-state u* described in Section HV A 21 

The simulation from which the snapshots in Figure [6] were extracted is very rapid, animating smoothly in real 
time at 30 frames/s even though N w 10'*, allowing immediate interactive model exploration. The raw calculation 
(no graphical animation) takes 0.008 s per time-step, benefitting greatly from the split operator method using FFT 
convolution. Performing this without the aid of the factorization (jlSp is about 100 times slower. Furthermore, we 
find that the additional effort needed to extract the mean-square displacement < > (t) := Jf^i^ ~ XQ)'^u{x,t)dx, a 
useful measure of spreading, is negligible. We remark that in order to model zero-flux boundary conditions while still 
using FFT convolution, a periodized kernel was used but the w array was zero-padded with a border of width 0{L), 
at negligible extra cost. 



V. LOCAL BEHAVIOR NEAR A SHARP HABITAT TRANSITION 



We now examine in more detail what takes place at the discontinuities in preference function w{x) that arise at the 
boundaries between different habitat types. As noted earlier, such discontinuities were not able to be treated in the 
derivation of the advection-diffusion limit in previous work [lO| . 

For simplicity, we consider a single discontinunity on a ID landscape, however we expect, and observe, similar 
behavior in the more biologically relevant case of multiple discontinuities arising at the edges of different habitat 
types on a 2D landscape. Consider a landscape in which there exists a single boundary between two habitat types 
located at a; = 0, resulting in the following discontinous preference function: 

wix) = I ^' ^ < [! (23) 

Given a jump kernel of width L it follows from ([7]) that z{x) is a smoothed (mollified) step- function with transition 
region width L. The analytic formula p3p then tells us that the steady-state pdf jumps by a ratio of wq precisely at 
the habitat transition; however, when viewed on length scales larger than L, the pdf jumps by a ratio Wq. For the ID 
exponential kernel p^ . the analytic expression for the steady-state follows from that of z via (O, and is 

-*(^) - { w[^'^e-^f^] , X J (2^) 

The spatial decay length is thus the same as that of the kernel </>. This is shown in Fig. [7^. 

How fast is this equilibrium reached? We demonstrate in Fig. [7)3 that the transition region's shape reaches its 
approximate equilibrium in only a single time-step (a slight change also occurs over the next few time-steps). Here 
we chose initial conditions which already matched the global pdf ratio of Wq, in order to study the equilibration in 
this local region alone. 



A. Effective boundary condition for advection-diffusion equation 

Armed with this understanding of the local behavior at a preference function discontinuity, we can incorporate this 
into the Fokker-Planck PDE model for the evolution of u{x, t) in the r ^ limit (in which case kernel width L must 
also go to zero at an appropriate rate). We remind the reader that the Fokker-Planck equation is 

|.^-|-[c(.).] + ^M(.).] (25) 

where c{x) and d{x), representing drift and diffusion rates, take on values given by the r-scaling of the second moment 
of the jump kernel (j) [2 [13, [ill 

We combine two observations: i) locally the steady-state in the vicinity of a discontinuity in w enforces a multi- 
plicative jump (the square of the w ratio) in u, and ii) in the r — > limit the Fokker-Planck equation evolves on 
much slower time-scales than the local equilibration in this vicinity (which happens in 0(r)). Thus we expect that. 
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for evolution on time-scales long relative to r the transition region is in local equilibrium, with the effective boundary 
condition 



for all t > 0, (26) 



where the subscripts — and + indicate limiting values on the left and right side of the discontinuity respectively. 
Similarly, by conservation of the flux J{x) = ■^[d{x)u] — c{x)u across the discontinuity, we must have that if d{x) is 
continuous and c{x) = then 

d d 

—u(x.,t)^—u(x+,t), foraUt>0, (27) 
ox ox 

at a step discontinuity in w. We may now interpret (j26p and (|27p as matching conditions for coupled advcction- 
diffusion equations on either side of the discontinuity. In this way we have a recipe to understand the diffusion limit 
even in the presence of discontinuous preference functions. We remark that our assumption of locally vanishing c 
corresponds to no gradient in w locally on either side (the case with general values of w' either side we postpone for 
future work). The above argument is non-rigorous, relying on reasoning based on separation of length- and time-scales. 
However it seems to be supported qualitatively by the evidence in Fig. [SjD-e (although at early times does not 
appear to hold accurately). We suggest that a more detailed analysis via matched asymptotic expansions should be 
carried out. 



VI. CONCLUSIONS 



We have analyzed the mathematical properties of a mechanistic resource selection (MERSA) model that captures 
the influence of spatially-localized habitat preference on the movement behavior of individuals and predicts their 
resulting patterns of space use. The model combines random foraging motion with a local sensitivity to habitat 
preference over a 'perceptual radius' L. Directed movement bias is generated in a similar manner to the angle-biased 
(von Mises) jump kernels used in [ll| . and also becomes equivalent in the small-L limit to continuous-time advection- 
diffusion in a 'confining potential' as in flj. Our analysis shows that the model has a desirable factorization (|15p which 
yields a simple closed-form formula for the steady state pdf. The effect of the compounded random movement 
decisions upon this steady state pdf is a novel geometry-dependent scaling with the preference function: linear when 
L is large (compared to local habitat features), but quadratic when L is small. 

These novel spatial effects are absent in conventional RSA, and have only been analysed previously in mechanistic 
home range models in the context of the advection-diffusion limit [101] in which the perceptual radius of individuals 
is small {L —> 0). Our analytic formula developed here allows this to be understood for the case of discontinuous 
preference functions and for any given perceptual radius L. 

Large gains in computational efficiency have been demonstrated throughout, including the case of discontinuous 
2D preference functions motivated by observations of spatially-varying prey availability across different habitat types. 
We believe such efficient forward models will be important tools, as inverse modeling, and fitting of multiple model 
parameters to observations, become more popular and time-consuming. 

Although we did not construct or study them numerically in this work, we expect that the benefits demonstrated 
here for a model that is linear in u (with prey- or spatially- dependent movement rates) will also apply in the analysis of 
non-/mear models such as those with density-dependent diffusion (e.g. Sec. 3 of [3), and scent-mediated interactions 
between multiple animal packs [l,[lll, [3. For example, chemotaxis could be included in the preference function w, in 
which case our analytic formula could be used to turn a coupled PDE system into coupled algebraic systems, a 
huge simplification. We also expect that by extending our preliminary operator analysis (Proposition [2]), the spectral 
properties, and hence equilibration rates of animal home range space use, may be deduced. 

Concise Matlab codes for computation of all figures in this work are freely available at 
http : //math . dartmouth . edu/~alib/moorcrof t/ 
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FIG. 4: Steady-state 2D pdfs u* for the exponential jump kernel (|2Up . applied to preference functions w derived from biomass 
data shown in Fig. |3l The three columns represent the choices a — 2, 15, 500 in the preference model (|22|l . At the top of each 
column is a surface plot of the preference function, followed below by surface plots of u" for three decreasing values of L. 
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FIG. 5: a) Piecewise linear model preference function wix), steady state pdf and predicted limiting case of steady-state 

pdf (dashed). b)-e) Snapshots of time evolution of u{x,t) under master equation Q in ID, at four times (indicated by 
t/r the number of iterations). The exponential jump pdf of (|19|l is used with L = 0.01. In e) u is very close to steady-state; 
the prediction (|lip is also shown (dashed). 
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FIG. 6: Density plots snapshots of the time evolution of u[x,t) under master equation ((2)1 in 2D, at four times (indicated by 
t/r the number of iterations). The exponential jump pdf of (|20p is used with L = 0.7, and preference function derived from 
biomass data in Fig. |3]via Ea. (|22|l with a = 500. In d) u is very close to steady-state, which is shown in Fig. |4]i. 
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FIG. 7: a) Detail of transition in steady-state space use occurring within distance L of a discontinuity in preference function 
of the form (|23p with wo — 2, for the ID exponential jump kernel of (|19|l . b) Time-evolution of u{x, t) in this same local region 
for 10 time-steps of the master equation ([2]), starting from an initial pdf ^0(2;) = w{x)'^. 



